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ABSTRACT 

A cellular automaton model of pulsar glitches is described, based on the superfluid 
vortex unpinning paradigm. Recent anal yses of pulsar glitch d ata suggest that glitches 
result from scale- invariant avalanches (Melato s et al.l 120071 ). which are consistent 
with a self-organized critical system (SOCS). A cellular automaton provides a 
computationally efficient means of modelling the collective behaviour of up to 10 16 
vortices in the pulsar interior, whilst ensuring that the dominant aspects of the 
microphysics are not lost. The automaton generates avalanche distributions that are 
qualitatively consistent with a SOCS and with glitch data. The probability density 
functions of glitch sizes and durations are power laws, and the probability density 
function of waiting times between successive glitches is Poissonian, consistent with 
statistically independent events. The output of the model depends on the physical and 
computational paramters used. The fitted power law exponents for the glitch sizes (a) 
and durations (b) decreases as the strength of the vortex pinning increases. Similarly 
the exponents increase as the fraction of vortices that are pinned decreases. For the 
physical and computational parameters considered, one finds —4.3 a ^ —2.0 and 
—5.5 ^ b ^ —2.2, and mean glitching rates in the range 0.0023 ^ A ^ 0.13 in units of 
inverse time. 
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1 INTRODUCTION 



It has recently been show n that glitches in ind ividual pulsars 
obey avalanche statistics (M elatos et al.1120071 ) . Such statisti- 
cal behaviour, in conjunction with the wide dynamical range 
of glitches, suggests that the underlyin g physics is that of a 
self-organized critical system (SOCS) |Pines fc Alparl 1 19851 : 



iBakll 19961 : |j nil 19981 : 1 Melatos et al l2007|). SOCS are char 
acterised by transitio ns between m etastable states via scale 
invariant avalanches !|jensenlll998h . The fractional increase 
in spin frequency spans seven decades across the entire glitch 
population (10" 11 ^ 5v/i> ^ 10 -4 ), and up to four decades 
in any individual pulsar. A natural way to model a scale 
invariant, avalanching system is to construct a cellular au- 
tomaton that is driven slowly to a threshold of instability 
|Field et al.lll995l ). Such model s have been explored in detail 
in the context of earthquakes jlSornette et al.lll99ll). g r anu- 
lar assemblies llBak et al.lll988l; IWiesenfeld et al.lll989l ; iBakl 
1 19961 : iFrette et al.l Il996l: IFruessner fc Jensen! 120031 ). solar 
flares (|Lu fc Hamiltonlll99ll ). and superconducting vortices 



Jjensenl Il99d: iField et all Il995l: [ Bassler fc PaczuskH Il998l : 
iLinder et al.ll2004 ). 

Theories of pulsar glitches centre around the mass 
unpin ning model first proposed by Anderson fc Itohl 
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Jl975l) and extended by many others (AlparetaL 
iPines fc Alpaill 19851 : lAlpar et all 19891 : iLink fc Cutlej|2002l ). 
The neutron superfluid in the stellar interior is threaded by 
many (~ 10 16 ) vortices, approximately one percent of which 
are pinned to the st ellar crust at grai n boundaries and/or 
nuclear lattice sites l|Alpar et al.| [T989). As the pulsar crust 
spins down electromagnetically, a lag builds up between the 
velocity of the pinned vortex lines (corotating with the crust) 
and the superfluid. When the transverse Magnus force (di- 
rectly proportional to the lag) surpasses a threshold value 
(equal to the strength of the pinning force), a catastrophic 
unpinning of vortices occurs, transferring angular momen- 
tum to the crust. In order for this mechanism to generate 
glitches on the scale observed, it requires up to 10 12 vortices 
to unpin simultaneously, exhibiting a high level of collective, 
non-local behaviour. 

The mass unpinning model spawned much activity 
in quantifying the micr ophysics that would lead to such 
macroscopic behaviour dPines fc Alpail Il985l: Ijonesl Il997l 
Il998al lbl: lPe Blasio fc Elgarovl Il999l : lElgarpy fc De B lasio 
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2001; Donati fc Pizzocherd 120031 ; lAvogadro et"ail \200H \. 

In particular, the strength of the pinning force has re- 
cently been discussed in detail using a self-consistent, 
semi-analytical mode l for the vortex-nucleus interaction 
l|Avogadro et alj 120071 ). yielding the unexpected result that 
the strongest pinning occurs in regions of the crust that 
are of high and low density, rat her than intermediate den- 
sity (|Donati fc Pizzocherd 20061 ). An interest has also been 
taken in the evolution and morphology of the neutron 
star crust, with particular attention paid to impurities 
and defects, which determine the density of pinning cen- 
tres available to the vorti c es, and their relative strengths 
l|De Blasio fc Lazzari|[l998h . lLink fc Cutler! (|2002h looked at 
the role that precession, if present, plays in the unpinning 
process, and discussed the dependence of the pinnning force 
on various microphysical parameters of the superfluid vor- 
tices and of the nuclear crustal lattice. 

Several other theories of pulsar glitches, similarly fo- 
cused on the microphysics of the superfluid in the stel- 
lar interior and crust, are currently under consideration. 
These include thermally driven glitches arising from sud- 
den change s in the frictional co upling within the neu- 
tron star jLink fc Epstein! Il996h. superfluid turbu lence 
i|Peralta et all 120051 . l2006bl lat IMelatos fc Peraltal 120071) . su- 
perfluid two-stream inst abilities ( Andersson et al.(~ 2004 



Mastrano fc Meratosll20o"5T). c r ust cracking llRudermarJ 19691 ; 
Bald 1 19961 ; lAlpar et ail 1 19961 iMiddleditch et all l200rj) and 



accretion (|Morlev fc Schmidtll 19961 ; I Morlevl 119961 ) . 
Regardless of the model one chooses to examine, there 
is a common need to synthesise the microphysics with 
the global, collective dynamics in order to accurately de- 
scribe the observational data. In order to produce a realis- 
tic model of the global dynamics, some subtle aspects of 
the microphysics need to be sacrificed. Cellular automa- 
ton models of complex systems aim to employ only the 
dominant aspects of the microphysics to generate replica- 
ble and effici ent automaton rules. Suc h an approach was 
attempted bv lMorlev fc Schm idt (1996) using a cellular au- 
tomaton platelet model of mass accretion onto the stellar 
crust. Although the microphysical details differ from those 
treated in this paper, the underlying philosophy of using 
a simple model to describe large scale behaviour is iden- 
tical to ours. We are similarly encouraged by the marked 
parallels between superfluid vortices and those that pop- 
ulate hard, type II superconductors. The agreement be- 
tween cellular automaton models of flux creep and vortex 
avalanc hes in superconductors and exper i mental data is ex- 
cellent llLinder et alll2004|; I Jensen! Il990l; iField et all h 9951 ; 
iBassler fc Paczuskil 1 19981 ; iNicodemi fc Jensen! l2001al bl), in 
particular with respect to the measurement of a power law 
over several decades in the avalanche (cf. glitch) size. 

A general, first principles theory of self-organized criti- 
cality does not yet exist for any SOCS, let alone for pulsar 
glitches. In this paper we look for empirical agreement be- 
tween the output of our cellular automaton and the gross 
qualitative features in pulsar glitch data. A detailed quanti- 
tative comparison between a first principles theory and data 
is not possible at this stage. 

In this paper we present a cellular automaton model 
of pulsar glitches based on the mass unpinning model. After 
describing the key features of a SOCS we revisit the observa- 
tional pulsar glitch data in S[2] Sj3] reviews the vortex unpin- 



ning paradigm and recent advances in the understanding of 
neutron star structure which are relevant to the paradigm, 
^describes in detail the cellular automaton model and jus- 
tifies physically the automaton rules. In S[5] the statistical 
behaviour of the model is explored with respect to driver 
parameters, internal physical parameters, and variations in 
the automaton rules. Further discussion of the model's va- 
lidity and a summary of our findings are contained in [j6] 



2 SELF-ORGANIZED CRITICAL SYSTEMS 
2.1 General properties 

A system is described as a SOCS based on two underlying 
behavioural patterns. Self- organized implies the ability to 
develop structures and patterns in the absence of manip- 
ulation by an external agent. Critical implies that, like in 
phase transitions near the critical temperature, a small lo- 
cal perturbat ion can propagate throughout the entire system 
i|Jensedll998l ). 

For a system to be in a self-organized critical (SOC) 
state it mus t first satisfy the fo llowing three conditions 
l|jensenlll998l ; IMelatos et aill2007l ): 

(1) It is composed of many discrete, mutually inter- 
acting elements, whose motions are dominated by local 
(e.g. nearest-neighbour) rather than global (e.g. mean-field) 
forces. 

(2) Each element moves when the local force exceeds a 
threshold. In this way stress accumulates sustainedly at cer- 
tain random locations (metastable stress reservoirs) while 
relaxing quickly elsewhere. 

(3) An external force drives the system slowly, in the 
sense that elements adjust to local forces rapidly compared 
to the driver time-scale. Combined with local thresholds, 
this ensures that the system evolves quasistatically through 
a history-dependent sequence of metastable states. 

If the above conditions are met, the fol l owing proper- 
ties are generically observed (|jensenl [l99sl ; IMelatos et ail 
EM): 



(4) Transitions from one metastable state to the next oc- 
cur via avalanches: spatially connected chains of local equi- 
libration events, in which one element relaxes and redis- 
tributes some local stress to its neighbors, which in turn 
can exceed their thresholds and relax (knock-on effect). 

(5) Avalanches have no preferred scale: they can involve 
a few (commonly) or all (rarely) of the elements in the sys- 
tem. Their sizes and lifetimes follow power law distributions, 
whose exponents are related. Numerical values of the expo- 
nents depend on the spatial dimensionality of the system 
the s ymmetry and strength of the lo cal forces dField et al 



1995), and the level of conserva tion (|Vespignani fc Zapperi 
1998 ; (Pruessner fc Jensenll2002l ). 



(6) Over the long term, the system tends to a critical 
state, which is stationary on average but fluctuates instan- 
taneously. 

A driving force is not essent ial to the operation of a 
socs. H rucssner fc Jensed (|2002l) assert that an external 
driving force compensates for a lack of conservation in the 
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microscopic dynamics of a SOCS and hence is superfluous if 
the system is conservative on microscopic scales. We claim 
that our model is locally non-conservative and hence a driv- 
ing force is essential to maintaining criticality (see 94. 7[) . 

Avalanches result from the relaxation of isolated 
metastable stress reservoirs. The storage capacity of these 
reservoirs is scale invariant. Each relaxation event is statis- 
tically independent from other relaxation events in size and 
when it occurs. Two or more simultaneous relaxation events 
at different locations cannot be distinguished by their global 
effect (e.g. spin up of a neutron star), which is the cumu- 
lative result of multiple relaxations. Likewise, relaxation of 
two or more smaller reservoirs causes an avalanche which 
is indistinguishable from that produced by the relaxation 
of a single large reservoir. Statistical independence implies 
that the waiting times between avalanches follow Poisson 
statistics, and this is confirmed by models based on cellular 
aut omata. 

lAlpar et al.l (|l995l . 1 19961 ) discussed the formation of de- 
pleted and capacitive regions of vortices in a neutron star 
interior, generated by steep gradients in the vortex pinning 
potential. Wherever the pinning potential is greater, the vor- 
tex density is also greater, resulting in a stronger local inter- 
vortex interaction (Magnus force). When the Magnus force 
exceeds the pinning threshold, vortices are redistributed lo- 
cally into neighbouring regions, relaxing the capacitive re- 
gion (stress reservoir). This qualitative picture bears the 
hallmark of a SOCS. 



2.2 Pulsar glitch statistics 

iMelatos et al.l (|2007l ) compiled all of the available glitch data 
for pulsars, analysing in detail the nine pulsars that have 
been observed to glitch at least five times. We consider here 
the subset that have glitched more than six times (N s ^ 6). 
Of the eight pulsars satisfying this criterion, six are consis- 
tent with a power law probability density function in the 
fractional glitch size Sv/u, where v is the spin frequency of 
the pulsar, implying a cumulative probability of the form 
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where (<5^/^) max is the maximum observed glitch size, 
and (8v / v)m\ -n is the mimimum observed glitch size. 
IMelatos et alj (|2007l ) found that a takes a different value 
for each pulsar; when the data from all nine pulsars with 
N s 6 are considered, Eq. ([1} with a single value of a is 
not a good description of the aggregate Sv/v distribution. 
Two pulsars also exhibit quasiperiodic behaviour, which is 
also a natural manifestation of avalanche dynamics; when 
the global forces (driven externally) overwhelm the nearest- 
neighbour interactions, quasiperiodici ty in the waiting time 
distribution is expected l|jensenlll99Sf ). 

For each of the six pulsars analysed here, we find all val- 
ues of the exponent a, for which the data are not inconsistent 
with the analytic distribution at the la confidence level ac- 
cording to the Kolmogorov-Smirnov (K-S) test. The range 
of a is presented in Table [1] Table [1] also gives the expo- 
nent that minimises the K-S statistic, along with the smal l- 
est glitch (Sv/u)mi n , for each pulsar l|Melatos et al.l 120071 ). 
However, it should be noted that all values of a within the 



la range (bounded by a m in and a max ) are equally likely; the 
exponent that minimises the K-S statistic is not a best fit. 
To appreciate why, suppose that the true underlying dis- 
tribution of 5v/v is known for some pulsar. Each time the 
pulsar glitches it samples this underlying 8u/v distribution. 
The glitch data for this pulsar represent one possible real- 
isation of the distribution containing N s glitches. There is 
no reason to attach greater significance to the a value that 
minimises the K-S statistic with respect to the particular 
real isation than to any o ther in the la range. 

IMelatos et al.l (|2007t ) also investigated the distribution 
of waiting times At, between successive glitches. Based on 
the statistical independence of each avalanche and the as- 
sumption that the system is driven at a constant rate, it 
is expected that the pro bability density function for At is 
Poissonian (Jensen 1998), implying a cumulative probability 
distribution of the form 

1 ^ exp(— AA£ m i„) — exp(— AAt) 
N g £-< exp(— AAt m in) - exp(-AAt max ) ' 



P(At) 



(2) 



The sum in © is taken over the minimum observable wait- 
ing time At m in, which is set by the gap between data spans 
in which a glitch is localised. At m i n is a function of the ob- 
serving schedule and hence of epoch; e ffectively, therefore, it 
takes a different value for each glitch l|Melatos et al.ll2007l ). 
At max is the total interval over which a pulsar is observed. 

Figure [T] shows cumulative histograms of waiting times 
At, for the same six pulsars listed in Table Q] plotted agaist 
XAt, where A is the Poisson rate for each pulsar in units of 
yr _1 . Once again, A is chosen to minimise the K-S statistic 
for definiteness, but any value in the range A m i n A ^ A max 
would do just as well at the la confidence level. The observa- 
tional data is overlaid with the ideal cumulative Poissonian 
1— e~ AAi for comparison, taking At mln — and At max — » oo; 
that is, we do not correct for dat a gaps and a fin i te tot al ob- 
servation for simplicity, unlike in Melatos et al.l (|2007l ) . The 
agreement with the universal Poissonian is good for these 
six objects, with 0.35 yr _1 ^ A ^ 5.6yr~ x . Observational 
confirmation that pulsar glitch waiting times obey Poisson 
statistics upholds the prediction that the glitches arise from 
the relaxation of metastable stress reservoirs isolated by re- 
laxed zones. 

The observed distribution of glitch durations cannot be 
discussed at this t ime, as very few glitch es have been re- 
solved temporally (|McCulloch et alj[l990h . 



3 VORTEX UNPINNING IN A PULSAR 

The mass unpinning paradigm of pulsar glitches assumes 
that the inner crust is a dense nuclear lattice permeated by 
a superfluid that is threaded by many quantised vortices. 
The vortices interact repulsively via the Magnus forc^B (per 
unit length of the vortex line) given by 



1 The Magnus force is in fact a fictitious force, similar in nature 
to the Coriolis or centrifugal forces. A free vortex, in the absence 
of external forces, moves with the local superfluid flow. In order 
for a vortex to move with respect to the local flow, a force must 
be applied such that the total force per unit length of vortex line, 
pn X (vl — v s ), is sufficient to sustain the motion of the line vl 
with respect to the local flow. 
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Figure 1. Cumulative distributions of glitch waiting times At, for six pulsars with N g 6. The horizontal scale XAt is different for 
each pulsar; A is the mean glitch rate that minimises the K-S statistic for that pulsar as listed in Table [T] The solid curve represents the 
ideal, Poissonian, cumulative distribution, P(At) = 1 — e _AAt . The waiting time distribution in an ideal SOCS follows this Poissonian. 



Table 1. Parameters of the glitch size and waiting time distributions for pulsars with N s ^ 6 l|Melatos et al.ll2007l l ■ a and A are chosen 
to minimise the K-S statistic. 
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where p s is the superfluid density, n is a vector pointing 
along the vortex, whose modulus \k\ = h/2m n is the circu- 
lation around a single vortex (m n is the neutron mass), vl 
is the vortex line velocity, and v s is the local bulk superfluid 
velocity. If allowed to move freely, vortices organize into a 
rectilinear array spaced evenly according to Feynman's rule 
(vortex area density n = 4itv/k). In keeping with condition 



(1) in §2.11 each vortex is considered a discrete element of a 
SOCS. For the purpose of practical modelling, given that we 
are dealing with ~ 10 16 vortices in a neutron star, we group 
vortices into bundles (in internal equilibrium) and regard 
the bundles as discrete elements. The velocity field gener- 
ated by each vortex is inversely proportional to the distance 
from the vortex line, so local interactions dominate global 
forces. 

Condition |(2)| in §2.11 exactly describes the pinning of 
vortices to the pulsar crust, whether the pinning centres 



are nuclear lattice sites, defect s [monovacancies or impuri- 
ties (|De Blasio fc Lazzarilll998h ]. or grain boundaries. Con- 



dition (2) does not rule out pinning of neutron and proton 
vortices in the neutron star's core. Vortices are unpinned by 
the bulk super fluid fl o w wh en the Magnus force exceeds the 
pinning force. iJonesI (|l997n suggested that vortex pinning 
is too weak to occur in the pulsar crust, and that strong 
magnetic interactions between neutron and proton vortices 
are a feasible alternative. In the eve nt tha t strong p inning 
is provided by monovacancies alone, IJonesI l|l998al lbl) calcu- 
lated the monovacancy concentration needed to provide a 
pinning force greater than the Magnus force, and claimed 



that, for p B > 3.4 x 10 1 
physically large. 



kgnT 



the concentration is un- 



in Cniis the 



The external driver fulfilling condition (3) 
accumulating Magnus force generated by the build up of 
rotational lag between the vortex lines and the superfluid 
flow due to the electromagnetic spin down of the star. The 



A 



cellular automaton model of pulsar glitches 5 



effect of this slow build up is two-fold. Firstly, the vortices 
tend to slowly move radially outward, as the system strives 
to maintain the Feynman vortex density (jAlpar et al.ll 19841 : 
|jones|[r998al ). Secondly, abrupt adjustments occur when lo- 
cally the pinning force threshold is exceeded. Once a vortex 
is unpinned, it readjusts rapidly, either by returning to the 
'soup' of unpinned vortices (remembering that only ~ 1% 
of vortices are pinned at any one time), or by repinning 
at a nearby lattice or defect site. Most mesoscopic vortex 
models assume that the inertia of vortex lines is negligible: 
once unpinned, vortices move immediately with the bulk 
super fluid flow, regar d less of the d i rectio n of the unpinning 
force |Donnellvlll9'9l"h . iField et al.l (|l995T l pointed out that, 
at least in the context of flux creep in superconductors, such 
an assumption leads to 'stick-slip' dynamics that are more 
likely to resemble those of an ideal sandpile (the quintessen- 
tial SOCS), because inertial effects do not dominate the ef- 
fect of the threshold. The time-scale of the readjustment is 
much shorter than the driving time-scale of the stellar spin 
down. 

Throughout this paper we assume that superflui d vor- 
tices ar e straight and parallel to the axis of rotation. Ijonesl 
l|l998bl ) argued that the high, but finite, tension in vortex 
lines justifies this assumption, as the displacements caused 
by the pinning centres are much smaller than the nuclear 
separation. Contrarily, it is the finite tension in vortices that 
allows them to be drawn into the potentia l well of the pin- 
ning c entres and pin at all (jjonesl Il998bh . iLink fc Cutlerl 
(2002) suggested pinning centres in the nuclear lattice ef- 
fectively span approximately thirty lattice sites, due to fi- 
nite tension, where the internuclear spacing is ~ 10~ 14 m. 
Relative to the coarse graining scale of a practical cellu- 
lar automaton (linear cell dimension up to ~ 2 x 10 3 m), 
this deviation from straight vortices is clearly insignificant. 
However, global hydrodynamic simulations suggest that the 
array breaks up into a vortex tangle above a critical ro- 
tational lag, driven by the Donnelly-Glaberson instability, 
with a n et polarisation parallel to t he rotation axis [see for 
example iPeralta et all (|2005l . l2006al lbl)] . 



4 CELLULAR AUTOMATON 

The cellular automaton model presented in this paper takes 
the basic physical features of the inner crust-superfluid sys- 
tem and interprets them as simple rules. Much of the mi- 
crophysics that has emerged from detailed studies of pin- 
ning strengths, crustal structure, and turbulent flow (see 
references in ^TJ) has not been taken into account directly 
when constructing the automaton rules, but it is relevant 
for understanding how the system can be 'scaled up' (renor- 
malised) to the necessary coarse grain. In terms of reproduc- 
ing the large-scale, collective dynamics of the SOCS, these 
simple cellular automaton rules have many advantages. In 
particular, the simple rules allow for a large range of system 
sizes and parameter space to be modelled for long enough 
to obtain stationary statistics, without prohibitive compu- 
tational cost. Table [2] summarizes the nomenclature used in 
this and subsequent sections. 



4.1 Grid 

The simulation grid contains JVg rid cells, covering a circular 
'star' of radius R, representing a projection onto the equa- 
torial plane. The grid points are arranged in a rectilinear 
(square) array in order to easily identify nearest neighbours 
and to ensure that all cells are equal in size. In reality, the 
equilibrium configuration of superfluid vortices is a triangu- 
lar Abrikosov lattice. For the purposes of this model, we as- 
sume that the two configurations are indistinguishable when 
vortices are bundled in large groups. Vortex positions are 
restricted to the coordinates of the centre of each grid cell. 
As mentioned in there is an inherent assumption that 
the vortices are straight and parallel to the axis of rotation, 
and that the vortex dynamics are independent of the vortex 
length. Indeed, the intervortex force is calculated as a force 
per unit length. 

The number of vortices that occupy the grid is calcu- 
lated using the Feynman relation 

j v s ■ dl = kN v , (4) 

where the integral is taken around a closed path (e.g. a circle 
of radius equal to the stellar radius), and N v is the total 
number of vortices enclosed by the path. In the limit of an 
infinite vortex distribution, the equilibrium distribution is a 
rectilinear array. Although we have a finite array of vortices, 
we assume that they form a rectilinear array, such that the 
total number of vortices is 

N V (R) = 4n 2 uR 2 /K . (5) 

Once we know how many vortices should occupy the 
system as a whole, we bundle them such that, on average, 
there is one bundle per cell. Partly, this is done to reduce 
the number of discrete elements and keep the computation 
practical. However, it is also done because uncertainties in 
pulsar timing limit the smallest glitch we can observe (corre- 
sponding to one bundle) . For any given pulsar we model the 
observed glitch range: the minimum glitch size G m i n (what 
constitutes a glitch is discussed in §4.5p occurs when one cell 
unpins, and the maximum Gmax occurs when all cells unpin. 
This restriction ensures that our model does not produce 
glitches outside the observati onally imposed range. Thi s is 
in contrast with the model of lMorlev fc Schmidtl (|l996t ) in 
which the minimum glitch size is set by timing noise and 
the maximum by the total available energy. Alternatively, 
we can have an average of B > 1 bundles per cell. In this 
scenario, the minimum glitch size is less than the minimum 
observed glitch size, and hence we can disregard all glitches 
resulting from the unpinning of fewer than B vortex bundles. 
We choose our grid size such that 

AW = ~ • (6) 

Pulsar timing data suggests, via (JBJ, a range of grid sizes 
varying from Agrid = 10 to a maximum of Agrid = 100 Q. 
The number of vortices per bundle is then 

AW n = N v / (SAW) ■ (7) 

2 ^grid ' s * ne total number of cells in the square grid. However, 
only cells that lie within the stellar radius are actually available. 
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Figure 2. Schematic of the cellular automaton. The left side of the figure is a cartoon of the grid, illustrating how the force imbalance 
on a cell is calculated. At the cell centred on the tip of the white arrow, the supcrfluid velocity field is calculated as the sum of three 
components: the circular velocity due to the the pinned vortices residing in the grey shaded circle; the circular velocity due to the 
unpinned vortices in the grey region; and the summed velocity field due to vortex bundles in the eight nearest-neighbour and next- 
nearest-neighbour cells (bundles are represented as black dots). The right side of the figure describes how the force imbalance due to 
the nearest-neighbour vortices, -FnNi is calculated. Each vortex bundle generates a velocity field according to a right-hand rule, with the 
circulation vector pointing out of the page, which falls off inversely with distance from the centre of the bundle. 



In an astrophysical context, however, we cannot assume 
that glitches smaller than the observable minimum do not 
take place, and with equal or greater frequency. If the system 
is in a self-organized critical state, we expect glitches to oc- 
cur at all scales, within the bounds of the system. To resolve 
smaller (Sv/v)m\n in the model, we must fractionate into 
smaller bundles. In order to make comparisons with the ob- 
servational data, we preclude a priori these smaller glitches 
from occurring in our model. If the system is truly scale in- 
variant, this effective window function (coarse-grain scale) 
should not interfere with its overall statistical behaviour. 
In the flux cr eep model for superconductors developed by 
Uensenl Jl99Ch . each cell holds at most one particle. In our 
model each ce ll holds as many bundles a s the avalanche his- 
tory dictates. Bass ler fc Pacz uski ( 1998) take a similar ap- 
proach, allowing many point pins in an extended cell. 

4.2 Initial conditions 

Initially the vortex bundles are laid out at randorrjf] such 
that the mean occupancy of a cell is B bundles. Whether 
or not it is realistic to have a vortex distribution that is 
inhomogoneous on the scale of the grid cells (~ 2 x 10 3 m) 
depends upon the distribution of the pinning centres. This 
point is disc ussed in more detail i n f6] From the glitch data 
analysed bv lMelatos et all (|2007f ) . we hypothesize that our 
system is in a SOC state at any given epoch, with an inho- 
mogeneous configuration of pinned vortices, and take this to 

3 Vortex bundles are placed at a randomly selected cell, one at a 
time, until all jV(, un bundles have been distributed. 



be the initial state of the automaton. We do not model di- 
rectly how this inhomogeneous state arises from an initially 
homogeneous state because the latter configuration is inac- 
cessible to observations; it is disrupted by avalanches almost 
immediately after the neutron star is born. Given that the 
random initial conditions of the model do not necessarily 
define a critical state of the system, we allow the simulation 
to complete many (~ 10 4 ) time steps before considering the 
output, to ensure that criticality has been established. 

4.3 Pinning threshold 

The strength with which vortices are pinned to the stel- 
lar crust, Jthreah, plays an essential role in the mass unpin- 
ning model of pulsar glitches. Recent calculations imply that 
pinning to nuclear lattice sites is strongest in regions of low 
(~ 4.9x 10 15 kg m~ 3 ) and high density (~ 1. 6x 10 17 kgm" 3 ) 
with a maximum energy of E p ~ 3 MeV l|Avogadro et al.l 
l2007h . We take the pinning force to be the same everywhere 
on the grid, with value -E P /£, where £ ~ 10~ 14 m is th e su- 
perfluid coherence length (jElgarOv fc De B lasio 200l|); i.e. 
-Fthresh ~ 5 x lO^Jm -1 . For the purposes of the model we 
convert this to a threshold on the velocity lag, Av p i n — 
(v L - v s ) pin ~ 10 4 ms _1 . 

Vortex pinning at lattice defects is also possible. 
iLink fc Cutlerl (|2002l ) argued that if the pinned (coupled) 
fluid makes up 1% of the total moment of inertia, then pin- 
ning must occur throughout the crust, suggesting that the 
crust should be treated as an amorphous solid with ran- 
dom pinning sites. The randomness ensures that the con- 
tributions to the force of attraction towards the n uclear 
sites on either side of the vortex line do not cancel (IJonesI 
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Il998al ). Another possibility is to adopt the random pin- 
ning potential used to mo del flux creep in superconductors 
ijBassler fc Paczuski|[l998h . 

Vortices unpin when the relative velocity of the super- 
fluid and the vortex lines is large enough to produce a Mag- 
nus force Fm greater than the pinning threshold Fthrosh- If 
the Magnus force does not exceed Fthresh, the pinning force 
adjusts to balance the Magnus force exactly. Given that the 
vortices are assumed to be inertialess, the excess force be- 
yond the pinning threshold is of no relevance, as once the 
vortices unpin they flow with the ambient superfluid: 



F ■ 

-t pin 



Fm if Fm < -Fthresh 

otherwise . 



(8) 



In our model, the superfluid velocity at a given point 
is the sum of three components: the contribution from all 
pinned vortices that lie on or within the star-centred circle 
passing through the point (these vortices are assumed to be 
in a regular Feynman array on average); the contribution 
from all of the unpinned vortices lying within the above 
circle (also assumed to be in a regular Feynman array; i.e. 
mean-field approximation); and the contribution from the 
eight nearest-neighbour cells on the grid, as depicted in Fig. 

El 



v s, pinned ,un pin nod + V s NN 



(9) 



The first two contributions, hereafter named the mean-field 
terms, are calculated using Eq. Q, assuming that the ve- 
locity field is exactly tangential to the closed path (a good 
approximation which becomes exact in the thermodynamic 
limit N v — * oo). 

By assuming that the vortices are arranged in a regular 
array, we can evaluate @ with tangential v s , without know- 
ing the location of each vortex (the validity of this approach 
is verified in §5.1[) . The Magnus force is then calculated from 
Eq. @, with vl equal to the veloicty of the stellar crust and 
v s ,tot given by Eq. ((9|. Initially, the mean Magnus force is 
very small. As the star spins down, however, the pinned con- 
tribution does not diminish, while the unpinned contribution 
does, and the mean Av increases. 

We emphasise that the automaton represents a locally 
averaged model where each cell covers many pinning sites. 
Most of these sites are unoccupied because, for pinning at 
lattice defects (cf. seismic faults), the total number of vor- 
tices is much less than the number of pinning sites. 



4.4 Driver 

In this paper we consider two separate driving forces. Over 
long time-scales (> 8v/v), the subject of §5.71 the Magnus 
force ramps up as the lag builds up between the pinned vor- 
tices and the unpinned superfluid. Over short time-scales, we 
include thermally activated unpinning of the vortex bundles 
in random cells|j In both scenarios (short and long time- 
scales), small variations in Av = v_l — v s due to the eight 



4 The contents of a randomly selected cell are deposited in the 
eight neighbouring cells. A more realistic implementation may be 
to select the thermally unpinned cell only from those cells that 
are already close to the pinning threshold, and therefore most 
likely to unpin thermally. 



nearest-neighbour cells trigger avalanches in the system. By 
treating the nearest-neighbour interactions with greater pre- 
cision than the longer range contributions, we are empha- 
sizing the property outlined in condition (1) of $2] 



4.5 Operational definition of glitches 

In the automaton model a glitch corresponds to the unpin- 
ning of vortex bundles. Its size is determined by the total 
number of bundles (and hence vortices) unpinned. The as- 
sociated fractional spin up of the pulsar is then given by 
the fraction of the total moment of inertia carried by the 
recently unpinned fluid. For example, if each vortex bun- 
dle represents the unpinning of fluid carrying 10 -7 of the 
moment of inertia of the system, then an avalanche of ten 
bundles results in a glitch of size 5v/v — 10~~ 6 . 

In practice, the angular momentum transferred to the 
crust by vortex unpinning depends on the product of the 
number of vortices unpinned (equivalently, the reduction in 
superfluid rotation rate) and the moment of inertia of the 
region through which they move before repinning, or until 
the avalanche stops, reaching a new quasi-equilibrium state. 
The reason for this is found in the Onsager-Feynman rela- 
tion Q. The superfluid velocity at some distance from the 
rotation axis is proportional to the number of vortices en- 
closed, so the further the unpinned vortices travel radially, 
the more the superfluid decelerates, and hence the more an- 
gular momentum is transferred to the crust. The definition 
of glitch size given in this paper is therefore an approxi- 
mation, valid only when three conditions are met: (i) all 
vortices travel the same radial distance during an iteration 
of the automaton; (ii) avalanches peter out after travelling 
one or two grid cells radially; and (iii) there are no steep 
radial gradients in vortex density (trivially satisfied by our 
uniform automaton, but not necessarily realistic). We show 
in £ )5.4l that these conditions are fulfilled by our automaton 
primarily because the unpinning activity is restricted to a 
narrow annulus. However, a more sophisticated automaton 
which incorporates radial gradients is needed to address this 
issue properly. 

To facilitate future comparisons between the model and 
data, we clarify that the automaton describes spin-up events 
where the jump in spin frequency is unresolved by current 
timing experiments; t he maximum rise tim e consistent with 
existing data is ~40 s (|Dodson et al.ll2002l ). Such glitches are 
accompanied by several recovery time scales, the longest of 
which may extend to the next glitch, together with a simul- 
taneous jump in spin-down rate. The automaton applies to 
this sort of unresolved spin-up event, although it does not 
seek to model the recovery physics. Within the database of 
observed glitches, there are some glitches that do not con- 
form to all facets of the above definition. In general, these 
nonstandard glitches are observed in pulsars that have only 
giitched one or two times. They are not described by the 
automaton in its current form, but they are usually large 
and hence relatively rare, so their impact on the statistics is 
muted. Another kind of nonstandard glitch, which we make 
no attempt to model, is the time-resolved secondary spin up 
('after shock') seen in the Crab 20 to 40 days after a standard 
glitch (IWonget al.ll200lh . 
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4.6 Automaton rules 

The rules of our cellular automaton encompass the key fea- 
tures of SOC outlined in discrete elements dominated 
by local forces, a threshold above which the elements move, 
and a slow driving force. The algorithm that governs the 
dynamics of the model can be summarised as follows: 

(1) The grid is initialised by randomly placing bundles of 
vortices on the cells that lie within the star, as in §4.11 and 

(2) The force imbalance Fm — Fthresh is calculated at each 
cell based on the scheme outlined in §4.31 Cells where the 
force imbalance is positive are flagged as supercritical. 

(3) Half the vortex bundles in a supercritical cell are repo- 
sitioned on the adjacent cell that lies in the path of the bulk 
superfluid velocity vector, which may lie on the diagonal, 
leaving the supercritical cell with half its original contents. 



Table 2. Default physical and computational parameters used to 
analyse the model. 



Steps (2) and (3) are performed simulataneously, so that no 
bias arises from the order in which the cells are considered. 

(4) The number of vortex bundles unpinned is recorded 
cumulatively. 



(4) are repeated until no cells are supercrit- 



(5) Steps' 
ical. 

(6) The total number of vortex bundles unpinned is 
recorded as the glitch size for the current big time step 
(where 'big' is defined below). 

(7) (a) Either the vortex bundles occupying a ran- 
domly selected cell are redistributed into the neighbouring 
cells; or 

(b) the spin frequency of the crust (and hence 
vl/(27tJ?)) is decremented by uAt ta , the number of un- 
pinned vortices is proprotionally reduced, and the number 
of pinned vortices remains unchanged. 



(8) Steps (1 



(7) are repeated for a predetermined num- 
ber of time steps, corresponding to either an arbitrary period 
of time [if (7) (a) is used] or the number of time steps mul- 
tiplied by Afts (where At ta is the length of each time step), 



if (7)(b) is used 



Throughout this paper we refer to the completion of steps 
|(l)D(4)| as a small time step, and steps [(T)|j(5)| as a big time 
step. Roughly speaking, a small time step lasts for the time 
it takes a vortex to be advected from one cell to its neighbour 
by the local superfluid flow. Of course, this varies somewhat 
with position. A big time step corresponds roughly to the 
mean waiting time between glitches, but please note that 
an avalanche does not always occur at every big time step 
given the above rules. The total duration of the simulation 
depends on whether the driver is thermally activated vor- 
tex creep (J^ At ts <C v/v), or electromagnetic spin down 
QTAi ts > v/v) (see §T3j. 



4.7 Vortex motion 

The rules for moving vortices to an adjacent cell following an 
unpinning event rely on the (incorrect) assumption that |v s | 
is the same in the vicinity of all supercritical cells, such that 
the time taken to travel from one pinning site to the next is 
the same throughout the grid. That is, the discrete nature 
of the automaton means that the small time step sets the 
characteristic inter-cell travel time. This approximation is 
justified in two ways in the model. First, as long as the time 



Parameter 


Symbol 


Default value 


Grid dimension 


Wgrid 


100 


Stellar radius 


R 


10 km 


Initial spin frequency 


"0 


30.25 Hz 


Current spin frequency 


V 


13.90 Hz 


Spin down rate 


V 


OHzs" 1 


Total vortices in simulation 


N v 


4tt 2 i/R 2 / k 


Vortices per bundle 


N bun 


^v/iV g 2 rid 


Bundle factor (avg. bundles per cell) 


B 


1 


Pinning threshold 


F ■ 
1 pin 


10 4 ms- 1 


Pinned fraction 


e 


0.01 



to move one cell is much shorter than the driver time-scale, 
it is reasonable to set the small time step using the slowest 
vortices. If, as we claim in §5.11 most of the avalanche activ- 
ity takes place in a thin annulus, this approach is accurate, 
because |v 8 | is roughly uniform over the annulus. Second, 
only half the vortices in a supercritical cell are redistributed 
to its neighbours. The travel time for vortices at the far side 
of the supercritical (departure) cell (relative to the destina- 
tion cell) is greater than the travel time for vortices start- 
ing near the border. Crudely, if vortices are evenly spaced 
throughout the departure cell, about half of them have suf- 
ficient time to reach the destination cell during a small time 
step. As a corollary, the model contains a dissipative mech- 
anism, which ensures that avalanches eventually terminate. 
This fe ature means that our system is non-conservative. Fol- 
lowing IPruessner fc Jensen! |2002j) , we require an external 
driving force to maintain the system in a critical state. 



5 STATISTICS 

5.1 General features 

Three quantities in particular are used to characterise our 
cellular automaton model of pulsar glitches: the distribution 
of glitch sizes Sv/u, durations t, and waiting times At. The 
glitch duration is defined as the number of consecutive time 
steps at which there exist supercritical cells. 

The hypothesis that the automaton results in avalanche 
dynamics, driven slowly by electromagnetic spin down or 
random thermal unpinning, is now tested. The model is al- 
lowed to run for 10 5 big time steps. In §5T2f-§5T6l we consider 
time periods over which v does not decrease significantly. In 
§5.71 we consider longer time periods over which spin down 
becomes important. We characterise the model by explor- 
ing its response to changes in the controlling parameters, 
both physical and computational. The physical parameters 
are the spin frequency v, the pinning threshold Av p i„, and 
the fraction of vortices that are pinned e. The compuational 
parameters are the numbers of bundles A r b un = BN slid and 
grid cells iVg rid . 

Fig. shows the raw output of the model in the form 
of a time series of glitch sizes for a standard set of parame- 
ters (defined Table Unless otherwise specified, these are 
the parameters used to generate the model output. This 
particular set of parameters is chosen to balance the need 
for sufficient avalanches to generate good statistics versus 



computational efficiency. The vertical axis gives the size in 
terms of the fraction of the total superfluid moment of in- 
ertia that unpins during the avalanche. The inset zooms in 
on a segment of the time series spanning 10 3 big time steps. 
Although the simulation runs for a total of 10 5 time steps 
(without spin down), the first 10 3 time steps are di scarded 
to allo w the system to establish stationary statistics. IJensenl 
(1998) asserted that fractal dynamics necessarily arise in 
a SOCS. Although we do not endeavour to quantify this 
mathematically, we are able to demonstrate, by considering 
time series such as the one shown in Fig. [3] that our system 
produces self-similar dynamics; qualitatively, the time series 
looks the same whether we plot 10 5 or 10 3 big time steps. 

The lower panel of Fig. [3] shows the probability den- 
sity functions of glitch sizes 8v/u, durations t, and wait- 
ing times At, plotted as histograms, for the standard set 
of parameters defined in the caption. Both 8vjv and t 
obey power law statistics, with probability density functions 
p(8vjv) oc (8v/v)~ a and p(t) oc t~ b respectively. This is a 
characteristic feature of a SOCS (see §2.ip . For the standard 
parameters in Fig. [3] we obtain a = 2.9 and b — 4.8. More 
generally, for the parameter ranges studied in this paper, we 
find 2.0 < a < 4.3 and 2.2 < b ^ 5.5 (see jO HH} . In a 
SOCS, a and b axe related by (|jensenlll998h 

& = l + (a-l)— , (10) 

73 

such that the size of an avalanche scales with its linear extent 
L and (8u/u) oc L 72 and the duration scales as t oc L 73 . The 
exponents in Fig.[3]imply 72/73 ~ 2.0, which is exactly what 
one expects for compact, two-dimensional avalanches with 
72 = 2 and 73 = 1 (see §5.T[) . Even more encouragingly, we 
find (b — l)/(o — 1) = 72/73 ~ 2.0 in many of our numerical 
experiments where criticality applies, even though a and b 
individually cover wide ranges. 

Based on the data analyzed in lMelatos et al. I i|2007f ). a.s 
well as the likeness between our model and other cellular 
automata that display self-organized critical behaviour, we 
expect the distribution of waiting times to be exponential, 
in keeping with a Poisson process (statistically independent 
avalanches). The bottom right panel in Fig. [3] confirms this 
expectation: the waiting-time probability density function is 
p(At) oc e~ AAt , with A = 6.8 x 10 -3 in units of an inverse 
big time step. To further verify the property of statistical 
independence, we calculate the linear Pearson correlation 
coefficient relating the size of an avalanche and the size of 
the avalanche immediately preceding it, as plotted in Fig. [4] 
We find that for 10 5 avalanches, the correlation coefficient is 
—0.012, implying that there is almost no correlation between 
an avalanche and its predecessor. 

Contrary to the general idea that a SOCS tends nat- 
urally to its critical state, we are forced to fine tune the 
parameters to elicit avalanches. The system is particularly 
sensitive to e and |Av p j n |. Similarly, the exponents of the 
best power law fits to the glitch size and duration distribu- 
tions (see §5.21 and §5.6|l depend on the same tuned param- 
eters. This is also observed to be the case experimentally 
when superconducting vortices r espond to a chang e in the 
ramp rate of the mag netic field ([Field et al.lll995h . In cel- 
lular automaton models of forest fires, the number of trees 
between two ignitions needs to be tuned according to sys- 
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Previous avalanche size (xlO 8 ) 

Figure 4. Next avalanche size versus previous avalanche size for 
9 X 10 4 avalanches, r is the linear Pearson correlation coefficient. 
The automaton shows no sign of 'memory' as \r\ <C 1. 



tern size, to avoid a cutoff or bump in the size distribution 
(IPruessner fc Jensen! [20021 ). 

Before analysing the results produced by the cellular 
automaton, we assess the accuracy of the mean-field approx- 
imation outlined in 34.31 Fig.[6]shows the velocity imbalance 
across the grid (as both a vector field and along a linear cut) 
generated using the mean-field approximation for a random 
distribution of vortex bundles (with a mean occupancy of 
one bundle per cell). The result is compared with an exact, 
'TV-body' algorithm, in which Av at any cell is computed by 
summing the v s produced at that cell by every bundle on 
the grid. The exact algorithm is not used elsewhere in this 
paper due to computational cost. Discrepancies between the 
mean-field and exact calculations are greater when a line is 
taken vertically through the velocity field, rather than di- 
agonally. This stems from the discrete, rectilinear nature of 
the grid. The two calculations are also in greater agreement 
further from the origin. By considering Av in the vicinity of 
the critical circle, where the mean-field contribution to the 
Magnus force approximately equals the pinning threshold 
(vertical dashed lines in right panel of Fig. , we see that 
the nearest-neighbour and gridding fluctuations are compa- 
rable. Yet the effect of gridding on the dynamics turns out 
to be negligible, because the output of the model is a power 
law in 8v jv. If the fluctuations were dominated by gridding 
effects, we would not expect power law statistics. 

The bold circle overlaid on the velocity field in Fig. [6] is 
called the critical circle. Well beyond this circle, the mean- 
field contribution consistently exceeds the pinning thresh- 
old, such that the nearest-neighbour variations are insuffi- 
cient to alter the supercritical state of the cell. Out here, 
vortex bundles never pin. Similarly, well within this radius, 
the nearest-neighbour variations are insufficient to alter the 
subcritical state of the cell. In here, vortex bundles never 
unpin. Near the circle, however, the nearest neighbour fluc- 
tuations are pertinent; indeed, they cause the avalanches we 
observe. The location of the 'active' annulus thereby defined 
depends on the physical and the computational parameters 
that define the model. 
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Figure 3. Top: Automaton output for 9 X 10 



time steps. The radius of the star is R = 10 km, the unpinning threshold 



[A v p!nl = 1.0 X 10 4 ms" 1 , the neutron superfiuid density is p B = 10 16 kgm - 3 , the initial spin frequency is Uq = 30.25 Hz, and the current 
spin frequency is v = 13.9 Hz. The avalanche size is given as the fraction of the total fluid moment of inertia that unpins during the 
avalanche. The fraction of pinned vortices is taken to be e = 0.01. The main panel displays a time series of avalanche sizes beginning 
after 10 4 time steps elapse, in order to ensure that a self-organized critical state has time to emerge from the initially random conditions. 
The inset zooms in on 5 X 10 3 time steps. Bottom: Probability density funtions of 5v/v {left), durations t (centre) and At (right) plotted 
as histograms on a log-log (left and centre) and log-linear (right) scale, using 10 2 logarithmically (left and centre) and linearly (right) 
spaced bins. The same parameters are used as for the top panel. The size and duration distributions are fitted with power laws of the 
form p(8u/u) = (&u/u)~ a and p(8v/v) = (&v/v)~~ h respectively. The waiting time distribution is fitted with an exponential of the form 



p(At) 



where A is in units of the reciprocal of one big time step. All fitted functions are plotted as dashed lines. 



5.2 Dynamic range 

In §4.11 we motivate the choice of grid size by referring to 
the observed dynamic range of pulsar glitches. In any indi- 
vidual pulsar the maximum dynamic range of 8v/v is 10 4 
(for PSR J1704-3015), and the minimum is 10 1 (for PSR 
J0537-6910EI)- Hence, from ©, we run our simulation on 
grids ranging from A gr id = 10 to Af gr id = 100 (the physi- 
cal size remains unchanged at R — 20 km). By construction, 
the range of avalanches outputted by the model should not 
exceed G ma x/G m m (i.e. 10 2 for the 10 x 10 grid, and 10 4 for 

5 This is one of two pulsars that glitch quasiperiodically; the 
other is Vela. 



JVgrid = 100 grid). Fig. [7] confirms this expectation. The top 
panels are histograms of the avalanches produced by sys- 
tems with A^id = 100, A?g r id = 60 and A grid = 10 over 10 s 
time steps. The dynamic range is greatest for the largest 
grid (10 -6,1 ^ 5v/v ^ 10~ 4 ' 2 ) and smallest for the smallest 
grid (10~ 4,2 ^ 8v/v ^ 10 -3 ' 2 ). It is not expected that these 
dynamic ranges span the full range G m in 8u/v ^ G ma x, 
as many of the grid cells lie inactive, well within the critical 
circle discussed in S I5.1I It is, however, encouraging that the 
larger grid sizes do in fact produce a larger dynamic range. 

Each of the histograms in Fig. [7]is overlaid with a power 
law best fit of the form p(8v/v) oc (8v/u)~ a , estimated us- 
ing a least-squares algorithm, with the fitted exponent in- 
dicated on each plot. By the central limit theorem, the un- 
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Figure 5. Fitted power law exponents for the distribution of glitch sizes (*) and durations (□). The error bars represent the lcr uncertainty 
returned by the least-squares fitting algorithm. Top left: 10 < Af grid < 100. Top right: 9.5 X 10 4 ^ |Av pin | < 10.4 X 10 4 ms~ 1 . Bottom 
left: 0.0092 < e < 0.010. Bottom right: 12.73 ^ v ^ 15.91 Hz. Automaton parameters: TV id = 100, B = 1, 10 5 big time steps. Physical 



parameters: p s = 10 kgm 



30.25 Hz, R = 20 km. 



Vector field Linear cuts 




Figure 6. Left: A vector plot of the velocity imbalance Av = vl — v s , generated using the mean-field approximation described in 
Eq. (0. The inset zooms in on the central 10 X 10 grid cells. The solid black circle shows the position of the critical radius, where the 
mean-field contribution to Av is approximately equal to Av p i n . Right: Profiles of Av along a straight-line cut through the grid, as a 
function of distance from the origin, as calculated by the mean-field (black) and exact (TV-body) (grey) algorithms. The solid and dashed 
curves correspond to the vertical and diagonal cuts drawn in the left panel. The two vertical dashed lines show the position of the critical 
circle. The horizontal dashed line marks Av D ; n . 
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Figure 7. Probability density functions p(5v/i/) of avalanche sizes Sv/v, plotted as histograms on a log-log scale, using 10 2 logarithmically 
spaced bins, for 10 5 big time steps. The histograms are fitted with a power law of the form p(Su/u) oc (8v/v)~ a {dashed lines). Top: 
Grids with N SI ^ = 100 (left), 60 (centre) and 10 (right). All panels have bundle factor B = 1. Middle: Bundle factors B = 1 (left), 
10 3 (centre) and 10 5 (right). All panels have iVg rid = 100. Bottom: B = 10 5 and Af grid = 100 (left), 60 (centre) and 20 (right). All 
panels have B = 10 5 . The vertical grey line denotes the crossover scale. Physical parameters: |Av p ; n | = 10 4 ms — 1 , p s = 10 16 kgm — 3 , 
v = 30.25 Hz, v = 13.9 Hz, e = 0.01, R = 20 km. 



certainty in each bin is inversely proportional to the square 
root of its contents. The fitted exponent is graphed as a 
function of ./Vg^d in Fig. [5] with error bars taken from the 
least-squares algorithm. There is a trend in both a and b to 
increase with increasing iV gr id, in keeping with the increas- 
ing dynamic range with N g ^. When the simulation is run 
for 2 x 10 5 time steps, the dynamic range does, as expected, 



increase; however, it still falls well short of G m ax/G m in. The 
number of bundles available to unpin is restricted to a small 
fraction (~ 15%) of BN gTid that lie in the vicinity of the 
critical circle. To better match the desired dynamic range, 
we should select the grid size such that G ma x/Gmin cells lie 
in the active zone near the critical circle. 

The glitch sizes generated by the model should be 
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Figure 8. Probability density functions p(t) of avalanche durations t, for bundle factors B = 1 (left), 10 3 (centre) and 10 5 (right), plotted 
as histograms on a log-log scale, using 10 2 logarithmically spaced bins. The histograms are fitted with a power law of the form p(t) oc t~ b 
(dashed lines). Automaton parameters: -/V grid = 100, 10 s big time steps. Physical parameters: |Av p j n | = 10 4 ms~ 1 , p s = 10 16 kgm — 3 , 
u = 30.25 Hz, v = 13.9 Hz, e = 0.01, R = 20 km. 
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Figure 9. Probability density functions p(8v/v) of avalanche sizes Su/u, plotted as histograms on a log-log scale, using 10 2 logarithmically 
spaced bins. The histograms are fitted with a power law of the form p(&v/v) oc (8v/v)~ a (dashed lines). Top: Current spin frequency 
v = 15.9 (left), 14.7 (middle) and 13.1 Hz (right). Bottom: Pinning threshold |Av pin | = 9.5 X 10 3 (left), 9.9 X 10 3 (middle) and 
10.4 X 10 3 ms~ 1 (right). Automaton parameters: N^-^ = 100, 10 5 big time steps, B = 100. Physical parameters: p s = 10 16 kgm~ 3 , 
u = 30.25 Hz, e = 0.01, R = 20 km. 
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Figure 10. Probability density functions p(At) of waiting times At, plotted as histograms on a log-linear scale, using 10 2 linearly spaced 
bins, for 10 s big time steps. The histograms are fitted with an exponential of the form p(At) oc e~ AAt (dashed lines). Top: Grid size 
JVgrid = 100 (left), 60 (centre) and 10 (right). Bundle factor B = 1. Bottom: Bundle factor B = 1 (left), 5 X 10 2 (centre) and 10 5 (right). 
Grid size 7V gri d = 100. Physical parameters: |Av pin | = 10 4 ins -1 , p B = 10 16 kgm -3 , u = 30.25 Hz, v = 13.9 Hz, e = 0.01, R = 20km. 



treated as a indi cative only. Evidentl y, several of the pul- 
sars analyzed bv lMelatos et al.l (|2007t ) exhibit glitches sev- 
eral orders of magnitude smaller than those produced by 
the automaton. It is easy to elicit such small glitches from 
the automaton by adjusting F p i n or v — vq so that the ac- 
tive circle lies closer to the rotation axis, thereby enclosing 
fewer vortices. We have not explored these possibilities here 
as there is too much freedom in the choice of parameters 
at present. (See §5.51 for further discussion.) If pinning is 
strongest in two or more annuli within the crust, the model 
predicts that glitches fall into two or more size brackets, de- 
pending on wh ere the annuli a re loc ated. It should be noted, 
however, that iMelatos et alj l|2007l ) found no evidence for 
a bimodal distribution when analyzing the latest statistics 
from individual pulsars, whose size distributions are consis- 
tent with unbroken power laws. 

5.3 Vortex bundling 

An underlying assumption of the model is that the scale on 
which the pinning centres and vortices are coarse grained 
does not alter the behaviour of the output qualitatively. In 
this section, we explore the impact of maintaining the grid 
size whilst varying the mean number of bundles per cell B. 



The simulation is designed such that the smallest (largest) 
avalanche occurs when one (every) bundle unpins. In a scale 
invariant system, avalanches should occur below and above 
the observational limits. For this reason we probe the effect 
of increasing the number of vortex bundles (equivalent to 
decreasing the number of vortices per bundle, for fixed v) 
on the distribution of avalanches. 

The avalanche size distributions for different values of 
the bundle factor B are shown in the middle panels of Fig. 
[7] The histograms are overlaid with a power law best fit. For 
B = 5 x 10 2 and B = 1 x 10 5 there is an obvious 'bump' in 
the glitch size distribution, whose peak occurs at the same 
value of 5v/v for different B [log 10 ((5^/i/)bum P = —4]. The 
bump is also present in the duration distribution (see Fig. 
[H), and there is a clear turnover in the slope of the waiting 
time distribution (see Fig. [ID) , which may be interpreted 
as a crossover from power law statistics to exponential de- 
cay. Importantly, crossover is evidence that the system is 
no longer in a SOC state, as there is a preferred scale for 
avalanches. 

In SOCS generally, there is a system-size-dependent 
crossover scale above which the power law distribution of 
glitch sizes gives way to exponential decay (see §2.I|I . In a 
genuine critical state, the crossover scale should be an in- 
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creasing function of t he system size (|jensedll998h . In the 
examples discussed by I Jensen] l|l998l ). the system size corre- 
sponds to the number of particles in the system (rice or sand 
grains, for example), not the physical extent of the system. 
That the bump does not appear for small B and becomes 
increasingly prominent for large B suggests a system size 
effect. On the other hand, the crossover scale, defined by 
the bump, does not increase with system size (i.e. B) in 
Fig. [T] Nor does system size explain why the bump is a 
local maximum in p(5v/v) inste ad of a turnover. In their 
experiments on superconductors iField et al.l 1 19951 ) saw a 
crossover to steeper decay in the avalanche size distribution 
for large values of the magnetic field ramp rate, but not a 
local maximum. If instead we interpret the system size to 
be the number of cells, not bundles, then the bump becomes 
more prominent for larger systems (e.g. with B = 1 x 10 5 ) 
as in Fig. [JJ Despite the bump not being present for smaller 
systems (7V gI .id 40), there is a value of Sv/v, for each iVg r id, 
above which p(8v/u) falls off faster than a power law, which 
is t he expected b ehaviour beyond the crossover scale defined 
by IJensenl l| 19981 ) (see Fig. [7J). Bumps like the ones in Fig. 
[JJ also arise in SOCS when one c hanges the morphology of 
the grains (e.g. rice versus sand ) (jFrette et al.lll996l ; I Jensen! 
ll99Sl ; |Pruessner fc Jensenll2003l ). 

The violation of scale invariance is at least partly due to 
a mismatch in how the pinning centres (cells) and vortices 
(bundles) are renormalized. The most scale invariant results, 
closest in spirit to a SOCS, are obtained for B = 1 (one bun- 
dle per cell, on average). By increasing B without changing 
-Wgrid, we are renormalising the vortices on a smaller coarse 
grain, whilst maintaining the coarsene ss of the grid, which 
is not a justifiable physical scenario. [Pruessner fc Jensen! 
(2002) claimed that inappropriately chosen parameters pro- 
duce a bump in the size distribution, which further suggests 
that B»l alters the fundamental statistical properties of 
the driven system. The implications for the microphysics of 
pinning in pulsars are explored in Sj6] 

In summary, bundling in the model introduce a pre- 
ferred scale into the system. (This scale greatly exceeds the 
vortex quantum which is too small to be resolved by a 
practical simulation.) When the average number of bun- 
dles per cell is kept near one, the system exhibits scale- 
invariant behaviour, as demonstrated by the power laws in 
glitch size. In this regime, the grid cells and bundles are 
well matched. If, however, the number of cells is significantly 
outnumbered by the vortex bundles, the scale invariance of 
the system is broken. In the absence of computational lim- 
itations, one should represent each pinning site and each 
vortex uniquely, eliminating the need for bundling. On the 
other hand, bundling confers a distinct advantage: it cir- 
cumvents the subtleties of vortex-vortex interactions and 
reconnections by renormalizing these effects into effective 
vortex- vortex forces on a coarsened grid. 

5.4 Pinning threshold 

The pinning threshold F p i n determines the minimum lag be- 
tween the superfluid and pinned vortices necessary to unpin 
vortex bundles. The pinning threshold can be reexpressed 
in terms of a minimum |Av p i n | between the pinned vortices 
(comoving with the crust) and the unpinned superfluid. 
As |Av p in| increases, the number of avalanches occur- 



ring within a fixed time (10 J big time steps) decreases, from 
33549 for [Av pln | = 9.5 x 10 3 ms _1 to 2227 for |Av pin | = 
10.4 x 10 3 ms~ 1 , as shown in Fig. [5] There are two rea- 
sons for this trend. First, greater imbalance is needed in 
the nearest neighbour cells to overcome the pinning force 
and trigger an avalanche. The greater is the required im- 
balance, the less probable it becomes. Second is the issue 
of fine tuning. For a given |Av p i n | there is a critical circle 
where the mean-field contribution to the velocity imbalance 
matches the pinning threshold almost exactly. Well outside 
the circle, every cell is supercritical and the vortices never 
pin. Near the circle, however, we find a mixture of sub- and 
supercritical cells because the nearest neighbour imbalance 
can act to reduce |Av| at a given cell, if it opposes the 
mean-field contribution. As |Av p i n | increases, the critical 
circle moves towards the origin, and hence the number of 
cells in its vicinity decreases. Fig. 1111 depicts images of the 
grid after 10 5 big time steps. The colours encode the value of 

Av| on each cell. Supercritical cells are white; all other cells 
are subcritical. The pinning threshold in the upper panel of 
Fig. [TT]is |Av p i n | = 10 4 ms~ 1 . We observe a smattering of 
supercritical cells near the critical circle. Inside the circle, 
all cells are subcritical. In the lower panel of Fig. 1111 with 
Av p i n = 9 x 10 3 ms _1 , we observe an annulus in which most 
cells are supercritical (where the vortices never pin). Just in- 
side the annulus, there is, like in the upper panel, a mixture 
of sub- and supercritical cells. It is in these mixed regions 
where we claim that self-organized criticality controls the 
unpinning dynamics. These results are evidence that the un- 
pinning activity is restricted to a narrow annulus, and hence 
unaffected by the use of a full circle as the simulation grid, 
as discussed in § 14.51 

A power law describes the size distribution well for all 
three values of |Av p i n | in Fig. [9] The best fit power law and 
its exponent are shown in each panel of Fig. [9] We find that 
a decreases as |Av p ; n | increases, however there is a slight 
upturn in both a and b for the largest values of Av p i n | con- 
sidered. For larger values of |Av p i n |, avalanches involving 
many cells (and hence many vortex bundles) are less likely 
to occur, both because the active annulus shrinks, and be- 
cause it is relatively unlikely that the nearest-neighbour im- 
balance sends a cell supercritical. Hence a decreases because 
the dynamic range decreases; the area under the histogram 
equals the number of avalanches. 

Fig- E] demonstrates that, for a finite, square vortex ar- 
ray, the profile of Av along a linear cut flattens out in the 
region where r ~ R, compared to smaller radii. Hence if 

Av p i n | approaches the maximum Av on the grid, more cells 
may participate in an avalanche, in comparison to when the 
critical circle lies further in. This may account for the upturn 
in a versus |Av p i n | seen in Fig. 

5.5 Pinned fraction 

The pinned fraction e impacts on the dynamics in two main 
ways: it changes the size of each vortex bundle, and alters 
the relative contribution to the velocity imbalance from each 
term in Eq. Q. Increasing the size of the vortex bundles in- 
creases the minimum avalanche size, whilst maintaining the 
maximum avalanche size. Larger bundles also enhance the 
role of the nearest neighbour cells. As more vortices pin and 
unpin, the first and third terms in Eq. ([9)l decrease with re- 
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Figure 11. Images of the velocity imbalance |Av| for |Av p ; n | = 10 4 ms 1 (upper) and 9x 10 3 ms 1 (lower). White cells are supercritical; 
the darker the cell, the smaller the Av. 



spect to the contribution from the unpinned vortices. Given 
that the unpinned vortices are homogeneously distributed, 
they are the foremost contributors to the mean-field critical 
circle discussed in §5.41 Decreasing their relative contribu- 
tion reduces the mean-field effect, emphasizing the nearest- 
neighbour, stochastic dynamics of the model. The lower left 
panel of Fig. [S] supports this argument: the la error in the 
fitted a and b values decreases with increasing e. The lower 
left panel of Fig. [5] also shows that both a and b increase 
with increasing e. For e ^ 0.0097 the a distribution flattens 
out, whereas the b distribution continues to increase. 

For e ^ 0.009 no avalanches are observed; the model 
is badly tuned. Moreover, for e > 0.010, the model slows 
down computationally, as a large number of cells regularly 
pin and unpin. This computational limitation arises because 
our algorithm does not excise regions in which vortices are 
permanently unpinned (see the lower panel of Fig. II 1 ft . in- 
stead calculating Av unneccessarily in these regions at each 
small time step. We intend to address this flaw in future 
work. 



5.6 Spin frequency 

Two values of the spin frequency enter the model: at birth 
(va), and today (y). As a first approximation, the total num- 
ber of pinned vortices is taken to be unchanged over the 
lifetime of the star and hence is dictated by vq as in Eq. ([5]). 
The pinned portion of the superfluid does not spin down 
with the stellar crust, whereas the unpinned superfluid does; 
its circulation is dictated by v. In the current model we 
hold v and vq constant (cf. §5. 7ft . with i/o — 30.25 Hz, and 
12.73 v < 15.91Hz. 

By decreasing v, we increase |Av| and hence the force 
imbalance at every cell. The impact is similar to decreasing 

Avpi n | because the critical circle shrinks. The two right- 
hand panels of Fig. [5] display a similar trend for increasing 

Avpi n | and v. The main difference is that v also determines 
the relative size of the first two terms in Eq. ([9)l , given that 
the number of pinned vortices remains constant while the 



number of unpinned vortices decreases with v. These effects 
increase the range of avalanche sizes as v decreases as shown 
in Fig. [9] For large values of v, the dynamic range goes to 
zero, as the critical circle has moved outside the star. For 
small v, a turnover emerges in the size distribution for small 
Sv/u, which is also present in the size distributions for small 
values of |Av p i n | in Fig. [5] This may happen because the 
critical circle lies well within the star, so there is a sizable 
region outside the critical circle where vortices are perma- 
nently unpinned. The automaton does not recognise that 
permanently unpinned vortices do not, in fact, participate 
in avalanches in real physical systems, so the output is bi- 
ased towards large avalanches. 

5.7 Long-term activity 

To model the long-term glitch behaviour of a pulsar, we must 
allow v to mimic the electromagnetic spin down of the stellar 
crust. Observed spin down rates differ by several orders of 
magnitude; in this paper we take v = —7 x 10~ 13 Hzs" 1 by 
way of example. Observations of pulsar glitches span a pe- 
riod of approximately forty years (~1970 to present), which 
equates to a total spin down of ~ 0.1 Hz. In our model, each 
big time step must therefore span ~ 10 6 s, in order to rep- 
resent a total observing time of forty years. Such a large 
time step determines the scale on which we can resolve indi- 
vidual avalanches. Ob servationally, glitches can be resolved 
down to a 40s window l|McCulloch et ail 199(1 iDodson et alj 
120021 ). 

Spin down causes the critical circle to shrink towards 
the axis of the star, because the number of pinned vortices is 
fixed, whilst the number of unpinned vortices decreases with 
v. The rate of spin down thus determines how quickly the 
glitch behavior changes as the critical circle moves through 
the strata of the star. Fig. [12] shows the size distribution 
(Su/u) for time steps At ts = 10, 10 5 and 5 x 10 5 s. Fig. [151 
shows the fitted power law exponents for the sizes and dura- 
tions, for 10 ^ Atts $S 5 x 10 5 . Within the uncertainty of the 
fitting algorithm, a and b are the same for all values of Atts 



A cellular automaton model of pulsar glitches 17 



10.4 



O 



At„ = 10 




-5.0 
•ogio^A) 



10.3 



At ls = 10 J 




-6.28 



-5.22 



-4.16 



10.55 



7.50 



4.45 



-6.28 



At ls = 10 7 




-5.12 
l°gio(<5fA) 



-3.96 



Figure 12. Probability density functions p(Su/u) of avalanche sizes &v/v, for time steps lasting Aits = 10 (left), 10 (centre) and 
5 X 10 5 (right), plotted as histograms on a log-log scale, using 10 2 logarithmically spaced bins. The histograms are fitted with a power 
law of the form p(8u/u) oc (8u/u)~ a (dashed lines). Automaton parameters: iVgjjd = 100, 10 5 big time steps. Physical parameters: 
|Av pin | = 10 4 ms-\ p s = 10 16 kgm- 3 , u = 30.25 Hz, v = 13.9 Hz, e = 0.01, R = 20km. 
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Figure 13. Fitted power law exponents for the distribution of 
glitch sizes (*) and durations (□), for time steps in the range 
10 $C Aits 55 5 X 10 5 . The error bars represent the la uncer- 
tainty returned by the least-squares fitting algorithm. Automa- 
ton parameters: N gI i^ = 100, 10 s big time steps. Physical pa- 
rameters: |Av pin | = 10 4 ms~ 1 , p s = 10 16 kgm~ 3 , u = 30.25 Hz, 
v = 13.9 Hz, e = 0.01, R = 20 km. 



tested. For the largest At ts , the time elapsed after 10 big 
time steps is 5 x 10 5 s, which equates to Au = —0.007 Hz. 
This moves the critical circle towards the origin by a negli- 
gible ~ 0.007%. 

The long term fate of a pulsar, with respect to its glitch 
behaviour, depends on the location of its pinning zones. If 
pinning can occur througout the pulsar, glitches continue to 
occur indefinitely, as the critical circle moves inwards. As 
time passes, however, the glitches become smaller, on av- 
erage, as the number of cells in the vicinity of the critical 
circle diminishes 0. If, as claimed by iDonati fc Pizzocherol 
(2003), pinning occurs at intermediate densities, no glitches 



occur once the critical circle shrinks inside the radius at 
which these densities occur. Alternatively, if pinning oc- 
curs at several strata [high and low densities, for example 
|Avogadro et al.ll2007l )]. glitches occur when the critical cir- 
cle lies in the vicinity of the pinning zones, which could mean 
that episodes of glitching are punctuated by 'quiet' periods. 

In the model s of post-glitch relaxation developed by 
lAlpar et al.l (|l996l ). the superfluid spins down when vortices 
move radially outward. The outward motion is driven by the 
distribution of capacitive and depleted zones; outward mo- 
tion is energetically favourable due to a radial bias in vortex- 
vortex scattering. In our automaton, spin down is an input; 
it does not arise self-consistently from vortex motion. By re- 
ducing the number of unpinned vortices with time (e.g. 5 15.70 
we attempt to recreate the scenario in which outward radial 
motion eventually makes vortices exit the pinning zones in 
the star. The automaton rules do inadvertently allow for 
outward radial motion thanks to the rectangular nature of 
the grid: azimuthal motion (i.e. in the direction of rigid body 
rotation) is never truly possible, as the vortices can only ever 
move horizontally, vertically, or diagonally. To incorporate 
outward motion and hence spin down self-consistently, the 
automaton rules should add a radial component to the bulk 
superfluid flow when computing vl- (This does not remove 
the artificial radial motion arising from the rectangular grid. 
In this context, an annular grid topology is preferable.) 



6 DISCUSSION 

In this paper we present a cellular automaton model for 
glitches in neutron stars, based on the vortex unpi nning 
paradigm l|Anderson fe Itohl Il975l ; lAlpar et al.l [l984). De- 
spite the gross idealisations in such a discrete model, the 
resulting statistical b ehaviour agrees w i th recent analyses 
of radio timing dat a dWang et al.1 |2000| ; IWong et all 1200 ll : 
iMelatos et al.ll2007l ). Of particular interest are the power 



6 Interestingly, there is no guarantee that glitches emanating statistical distribution (e.g. value of a). If so, we do not expect a 

from different critical circles are, in fact, drawn from the same pure power law in 5v/v over long times. 
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law distributions of glitch sizes and durations produced by 
the model, providing further evidence that glitches are a 
scale invariant phenomenon. For this reason, as observations 
improve in resolution and lengthen in duration, we expect 
to discover smaller and larger glitches than have been seen 
so far. Moreover, the exponential waiting-time distribution 
generated by the model supports the idea that the each vor- 
tex avalanche is a statistically independent event. The values 
of the power law exponents and Poisson glitching rates de- 
pend on the particular parameter set. Our fits return size 
and duration exponents in the ranges 2.0 ?J a ^ 4.3 and 
2.2 ^ b ^ 5.5 respectively, and mean glitching rates in the 
range 0.0023 ^ A ^ 0.13. Although we do not have a first- 
principles theory against which to compare our results, we 
claim that our cellular automaton is consistent with a SOCS. 

The operation of the cellular automaton relies heavily 
on the a ssumption that vor tices are inhomogeneously dis- 
tributed. lAlpar et al.l (|l996l ) suggested that crust cracking 
creates capacitive regions where many vortices pin simulta- 
neously. These regions do not permit a vortex current, so 
pinned vortices in these regions do not participate in out- 
ward vortex creep. Crust cracking is a nonuniform process, 
suggesting that the regions themselves are inhomogeneously 
distributed. Our model conforms qualita tively with this pic - 
ture. Of particular interest is the claim bv lAlpar et ail (|l995l ) 
that the neutron star crust is divided coarsely into depleted 
and capacitive regions. If verified, this claim suggests that 
the level of coarse graining necessary to render our cellu- 
lar automaton computationally tractable matches the level 
of inhomogeneity that is realistically presesnt in a neutron 
star. 

Of course, the opposite viewpoint is also viable: the 
crust forms a large, defect-free crystal, allowing the vor- 
tices to occupy a contiguous sequence of pinning positions 
and preventing the ty pe of inhomo geneity represented by 
our initial conditions (| Jones! Il998al ). Importantly, such a 
scenario challenges any model based on collective unpin- 
ning of numerous vortice s , not just our cellular automa- 
ton. Blasio fc Lazzaril (|l998h concluded that the crust 
is relatively free of microcrystalline structures and vacan- 
cies, with impurities, or point defects, the most likely dis- 
ruptions to the crystal lattice (faults and dislocations were 
not considered). Alternatively, if pinning occurs mostly at 
grain boundaries and c racks, as originally postulated by 
lAnderson fc Itohl j 19751 ). then inhomogeneity on macro- 
scopic scales is possible. 

Irrespective of the density and distribution of pinning 
sites in reality, the model assumes for simplicity that vortex 
pinning occurs throughout the pulsar. Several recent stud- 
ies have shown that pinning forces sufficient to oppose the 
Magnus force develop only in specific zones of the crust, e.g. 
regions of high and lo w, rather than intermediate, denisty 
|Avogadro et al.l 120071 ). Restricting the simulation to these 
optimal pinning zones obviates the need to allow for a range 
of pinning strengths if the zones are thin layers. A model in 
which pinning does not occur everywhere in the star leads to 
a different computational topology, e.g. an annulus, or sev- 
eral nested annuli. However, we choose to be conservative 
in this introductory paper by not imposing a grid topology 
a priori, in order to see if a preferred topology emerges by 
itself. In fact, it does; we find that a narrow active annulus, 
situated at the critical radius, participates in the avalanches. 



However, it is important to experiment with other topolo- 
gies in future work, because the active annulus we find in 
this paper may still be an artifact. If the active region lies 
near the rotation axis, an automaton executed in an annulus 
is free of the shortcomings introduced by the 'permanently 
unpinned' regions mentioned in §5.51 To use the full circular 
grid in a meaningful way, we must include pinning of proton 
and neutron vortices in the core of the neutron star, together 
with pinning of neutron vortices in the inn er crust. In this 
case, the additional physics of entrainment llLink fc Epstein! 
1 19961; iRuderman et all 1 19981 ; ISedrakian fc Sedrakianl Il995l ; 
lAndersson et al.ll2004l ) and induced electric currents implies 
a more detailed set of automaton rules, outside the scope 
of the model presented here. It sho uld also be noted tha t if 
the superfluid vortices are tangled l)Peralta et al.l [2006bl lah . 
rather than forming a regular Abrikosov array, these two- 
dimensional annular and circular topologies are clearly too 
restrictive. 

The locations of pinning zones also determine the long- 
term fate of a pulsar as it spins down. If pinning occurs 
everywhere, spin down increases the proportion of vortices 
with |Av| > |Av p in| which are always supercritical (outside 
the critical circle) and do not participate in avalanche dy- 
namics as they never pin. If pinning is restricted to one or 
more annuli, the glitch behaviour of the pulsar changes sig- 
nificantly when the critical circle moves into and out of the 
pinning zones. 

Clearly, our cellular automaton does not include the re- 
sponse of the pulsar to glitches (ie. the spin up of the pulsar 
crust and the subsequent 'relaxation'). Nor do we account for 
the finite time-scale on which the unpinned superfluid cou- 
ples to the crust, govern ed by the Ha l l-Vinen-Bekar evich- 
Khalatnikov equations (|Peralta et al.l 120051 . l2006bl ). We 
make the approximation that this time-scale is considerably 
less than the time-scale of our big time steps, motivated 
by glitch timing data (where the post-glitch recovery phase 
typically ends well before the next glitch; cf. Vela). 

To elicit avalanches from our model, we require fine tun- 
ing in both the physical and computational parameters, such 
that |Av max | ~ |Av p in|. This condition ensures that there 
are enough vortex bundles that switch between becoming 
sub- and supercritical as time passes. We achieve this by 
choosing v, e and Av p i n to place the critical circle near the 
surface of the star. We emphasize that for a larger (smaller) 
star, we would resize the critical circle proportinally, for 
computational rather than physical reasons. Fine tuning is 
also required in the ratio of vortex bundles to grid cells B. 
For B ^> 1, the automaton output is no longer consistent 
with a SOCS. 

In conclusion, we present an empirical cellular automa- 
ton model of pulsar glitches based on the mass vortex unpin- 
ning paradigm. We find that for certain physical and compu- 
tational parameters the model produces dynamics that are 
consistent with a SOCS and with radio timing data from 
pulsars. There exists no general, first-principles theory of 
SOC , let alone of pulsar glitches, so many of our results are 
empirical. In particular, ther is no way known at present to 
predict theoretically the size and duration distribution ex- 
ponents, and the mean rate of the waiting-time distribution. 
We do demonstrate that the basic physical principles gov- 
erning inter- vortex interactions can produce the type of col- 
lective behaviour necessary to explain pulsar glitches within 
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the mass unpinning paradigm, especially the puzzle of how 
so many vortices can unpin in sympathy during a glitch, and 
why their number varies so much from glitch to glitch. 

We thank Carlos Peralta for sharing his up-to-date 
glitch catalogue, and Stuart Wyithe for advice on statisti- 
cal methods. LW acknowledges the support of an Australian 
Postgraduate Award. 
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